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1  Summary 

Conventional  narrowband  FM  demodulation  approaches  assume  that  the  signal  of  interest  is  nar¬ 
rowband,  i.e.,  with  bandpass  spectral  content  confined  to  a  bandwidth  that  is  much  smaller  in 
comparison  to  the  carrier  frequency.  However  if  the  signal  of  interest  is  a  wideband  FM  signal,  the 
underlying  assumptions  are  violated  and  these  approaches  incur  significant  demodulation  error. 

In  this  report,  we  summarize  the  proposed  wideband  AM-FM  demodulation  approach  that 
employs  multirate  frequency  transformations  to  convert  the  wideband  FM  signal  to  a  narrowband 
signal  so  that  conventional  demodulation  algorithms  can  be  applied.  The  information  signals  asso¬ 
ciated  with  the  wideband  input  are  obtained  via  the  inverse  multirate  frequency  transformation. 

We  then  outline  an  approach  that  employs  multirate  Noble  identities  to  extend  the  proposed  ap¬ 
proach  to  larger  wideband  to  narrowband  conversion  factors  and  more  practical  implementations. 
We  further  detail  its  application  to  the  associated  problems  of  wideband  FM  image  demodulation, 
wideband  speech  formant  demodulation,  and  wideband  CPM  demodulation.  These  applications 
demonstrate  the  efficacy  of  the  MFT  based  wideband  demodulation  approach. 


2  Introduction 

Existing  conventional  narrowband  FM  demodulation  approaches  such  as  the  Hilbert  transform 
approach  [1]  or  the  Teager- Kaiser  operator  based  energy  separation  algorithm  [2]  assume  that 
the  information  bandwidth  of  the  communications  is  much  smaller  in  comparison  to  the  carrier 
frequency  of  the  signal.  When  this  condition  is  not  met,  specifically  when  the  signal  of  interest  is 
wideband,  they  incur  significant  demodulation  error.  Such  wideband  FM  signals  are  encountered 
in  satellite  communications  [3]  or  in  the  form  of  hypemasal  signals  encountered  in  people  with 
cleft  lip  or  palate  situations. 

In  recent  work  [4],  multirate  frequency  transformations  (MFT),  were  proposed  as  a  mechanism 
for  converting  wideband  FM  signals  into  narrowband  signals  so  that  conventional  demodulation 
approached  may  be  applied.  The  estimated  information  signals  are  then  inverse  transformed  to 
obtain  the  information  signals  associated  with  the  original  wideband  input. 

In  this  report,  we  outline  the  extension  of  this  approach  to  larger  wideband  to  narrowband 
conversion  factors,  thereby  enabling  more  practical  implementations  of  the  MFT  approach.  We 
further  detail  the  application  of  the  MFT  based  wideband  FM  demodulation  approach  to  the  asso¬ 
ciated  problems  of:  (a)  wideband  image  demodulation,  (b)  wideband  nasal  formant  demodulation, 
and  (c)  wideband  continuous  phase  modulation  (CPM)  demodulation  for  satellite  communications. 

2.1  AM-FM  Signal  Model 

Monocomponent  amplitude-modulation  frequency-modulation  (AM-FM)  signals  are  time-varying 
sinusoids  of  the  form: 


s(t) 


a(t )  cos  (  /  uji  (r)  dr  +  9\  )  , 


(1) 
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where  instantaneous  amplitude  (IA)  is  denoted  by  a(t)  and  the  instantaneous  frequency  (IF)  is 
given  by 

(f)  c  "l"  WmQi(t')  •  (2) 

Note  that  wc  is  the  carrier  (or  mean)  frequency  and  (pit)  is  the  normalized  baseband  modulated 
signal. 

Specific  for  sinusoidal  FM,  where  a(t)  remains  a  constant  A,  and  qt(t)  becomes  a  sinusoid,  the 
IF  can  be  further  expressed  as 


0Ji  {t)  =  c oc  +  Ujm  cos (ojft  +  6 2).  (3) 

Sinusoidal  FM  signals  can  be  expressed  via  Bessel  functions  as: 

+OO 

,(t)  =  A  J2  Jn  {P)  cos  (t oct  +  numt ) ,  (4) 

n=— 00 

where  Jn  is  the  nth  order  cylindrical  Bessel  function  of  the  first  kind.  The  modulation  index  of  the 
sinusoidal  FM  is  defined  as  the  ratio  [3  =  oom/uf  and  the  associated  carson’s  bandwidth  is  given 
by 

B  =  2(p+l)uf.  (5) 

If  1,  then  it  corresponds  to  the  wideband  FM  according  to  the  literature  of  FM  communi¬ 
cation  systems.  In  addition,  the  carrier-to-information-bandwidth  ratio  (CR/IB)  and  the  carrier-to- 
frequency  deviation  ratio  (CR/FD)  are  defined  respectively  as: 

CR  tuc  CR  coc 

T^T  =  — ,  -ER3  =  — ■  (6) 

In  l Of  r  JJ  um 

2.2  Conventional  Demodulation  Approaches 

In  this  section,  we  briefly  overview  three  conventional  demodulation  approaches.  Note  that  the 
adaptive  linear  predictive  IF  tracking  are  intended  for  FM  signals  only,  whereas  the  analytic  signal 
(AS)  approach  based  on  the  Hilbert  transform  (HT)  and  the  energy  separation  algorithm  (ESA) 
are  intended  for  AM-FM  signals. 

2.2.1  Adaptive  Linear  Predictive  IF  Tracking 

For  a  pure  FM  signal,  the  IF  demodulation  approach  employing  adaptive  filters  was  proposed  in 
prior  work  [4].  According  to  the  Wiener-Hopf  equations  [5],  the  optimal  coefficients  of  a  linear 
predictor  are  given  by 

=  R-^dx,  (V) 

where  wopt  denotes  the  optimal  tap  weight  vector,  R,r  denotes  the  input  correlation  matrix  and 
rdx  denotes  the  cross-correlation  between  input  vector  and  desired  signal.  As  summarized  in  [5], 
the  prediction  error  filter  is  given  by: 


E(z)  =  i-J2s7‘n, 

1=1 


(8) 
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where  x  are  the  coefficients  of  the  corresponding  optimal  predictor.  Conditioned  on  the 

small  prediction  error  assumption  and  the  further  assumption  that  the  message  signal  remains  es¬ 
sentially  invariant  over  the  sampling  range  of  the  linear  prediction  filter,  we  end  up  with  the  ap¬ 
proximation  in  [6]  through  (8)  given  by 

L 

^2gi(k)exp{-jl[wc  +  m(k)]}  ~  1,  (9) 

i=i 

where  gi(k)  is  the  weight  coefficient  of  tap  l  at  time  index  k  and  m(k)  is  the  sample  of  the  message 
signal  at  time  index  k.  Then  the  IF  of  the  signal  of  interest  can  be  estimated  by  executing  the 
following  steps:  1)  Compute  the  coefficients  of  the  prediction  error  filter;  2)  Obtain  the  roots  of  the 
coefficient  polynomial;  3)  Calculate  the  phase  argument  of  the  complex  conjugate  pole  location  of 
the  corresponding  roots. 

In  the  prior  work,  adaptive  algorithms  such  as  AS-LMS  (Analytic  Signal  Least  Mean  Square) 
and  AF-RLS  (Root  Least  Squares)  have  been  incorporated  into  the  structure  of  a  linear  predictor 
and  compared  with  each  other  based  on  the  demodulation  error.  However,  for  both  algorithms,  the 
step-size  or  the  forgetting  factor  needs  to  be  truncated  to  remain  within  certain  range.  Here  we 
choose  the  generalized  normalized  gradient  descent  (GNGD)  [7]  for  update  of  the  linear  predictor 
coefficients,  which  avoids  truncation  of  the  adaptively  adjusted  step-size.  The  algorithm  for  this 
GNGD  linear  predictive  filter  is  summarized  via 


e  (k)  =  x  (k  +  1)  —  ^  gi  (k)  x  (k  —  L  +  1) , 


i=i 


f3  ( k )  =  /3(k  —  1)  —  pa 
V  (k)  = 


e(k)e(k  —  l)xT  (k)x(k  —  1) 

W(m+mT2  ' 

a 


Mkni+my 

g  (k  +  1)  =  g  (k)  +  ri(k)e(k)x(k), 


(10) 

(11) 

(12) 

(13) 


where  x(/c)  and  g (k)  denote  the  vectors  of  input  and  tap  weights  at  time  index  k  respectively,  a 
is  the  step-size  parameter,  and  p  is  the  offset  learning  rate  parameter.  The  merit  of  the  GNGD 
algorithm  lies  in  the  adaptation  of  its  learning  rate  provides  compensation  for  the  assumptions  in 
the  derivation  of  normalized  least  mean  squares  (NLMS).  Due  to  its  robustness  and  stability,  the 
GNGD  is  well-suited  for  narrowband  nonstationary  signal  environments. 


2.2.2  Analytic  Signal  and  Hilbert  Transform 

For  a  one-dimensional  real- valued  signal  s{t)  :  R  — >■  M,  we  associate  with  it  an  complex- valued 
signal  z(t)  =  s(t )  +  jq(t),  where  the  imaginary  part  q(t)  is  defined  via 

q(t)  =  n  [s(t)]  =  s(t)  *\  =  [  —  (14) 

Kt  JR  tt4 

Note  that  q(t)  is  the  Hilbert  transform  of  s(t),  and  the  complex-valued  signal  z(t)  is  called  the 
analytic  signal.  Assume  that  s(t)  is  a  real-valued  AM-FM  signal  given  by 

s(t)  =  a(t)  cos  (c uct  +  4>{t)) , 
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(15) 


where  a(t)  is  the  IA,  t oc  is  the  carrier  frequency  and  0(f)  is  the  phase  for  the  IF  0'(f).  Note  that  the 
IF  defined  here  does  not  involve  the  carrier  frequency  c oc.  According  to  the  basic  property  of  the 
Hilbert  transform  and  Bedrosian’s  theorem,  q(t )  is  the  approximation  to  the  product  of  the  IA  and 
the  quadrature  of  the  FM  part  as: 

q(t )  =  H  [s(f)j  «  a(t)  sin  (t oct  +  0(f)) .  (16) 

However,  this  approximation  is  valid  only  if  the  following  conditions  hold: 

1.  a(t)  is  a  narrowband  lowpass  signal  that  varies  slowly  with  time, 

2.  carrier  frequency  c oc  is  sufficiently  large  such  that 

A  >  0'(t).  (17) 

As  a  result,  the  corresponding  analytic  signal  is  of  the  form: 

z{t )  =  s(f)  +  jq(t)  w  a(t)  exp  \j  (t oct  +  0  (f))] .  (18) 

Hence  the  IA  a(t)  and  the  IF  0'(f)  are  estimated  respectively  by 

a(t)  «  ||^(f)||,  (19) 

(20) 
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2.2.3  Energy  Separation  Algorithm 


The  energy  separation  algorithm  (ESA)  as  summarized  in  [2],  based  on  the  Teager- Kaiser  energy 
operator,  is  widely  used  for  monocomponent  AM-FM  demodulation.  For  a  continous-time  signal 
x(t),  the  nonlinear  Teager- Kaiser  energy  operator  in  the  continuous  case  is  given  by 

^c[x(t)}  =  x2(t)  -  x(t)x(t),  (21) 


where  x(t)  and  x(t)  denote  the  first  and  second  derivatives  of  x(t)  respectively. 

It  was  first  introduced  by  Kaiser  and  was  applied  to  track  the  energy  of  harmonic  oscillation. 
For  instance,  the  energy  of  the  oscillatory  signal  x(t)  =  Acos(coct  +  9)  can  be  tracked  by  the 
operator  \PC  via 

tyc[Acos(ujct  +  6)]  =  (Aloc)2.  (22) 

The  energy  operator  can  be  further  used  to  analyze  the  oscillation  of  signals  with  time-varying 
amplitude  and  frequency.  By  applying  the  energy  operators  \l/c  on  the  AM-FM  signal,  the  IA  a(t  ') 
and  IF  Ui(t)  (excluding  the  carrier  frequency  uc)  satisfy  the  following  relation: 

[a(t)  (uc  +  Ui(t))]2 .  (23) 


a(t)  cos  (  ujct  +  /  Ui(r)d,T  +  6 


As  a  result,  the  IA  a(t)  and  the  IF  cu;(t)  of  an  AM-FM  signal  x(t)  can  be  estimated  via  the  contin¬ 
uous  energy  separation  algorithm  (CESA)  summarized  by 


\/WW\ 


a(t)\. 


(24) 


Wi(t), 


(25) 


where  we  assume  that  the  IA  a(t)  and  the  IF  uj,  ( t )  do  not  vary  too  fast  or  too  greatly  in  value 
compared  to  the  carrier  frequency  c uc. 
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2.3  Narrowband  Constraint 


Conventional  demodulation  techniques  perform  well  only  under  the  narrowband  constraint.  In 
general,  the  IA  is  required  to  be  slowly-varying  and  the  IF  cannot  change  too  fast  or  too  greatly 
in  value  compared  to  the  carrier  frequency.  To  quantify  the  narrowband  constraint,  assume  the  IA 
a(t)  and  the  IF  uj,  (t)  are  both  bandlimited  with  highest  frequencies  ua  and  c Of  respectively.  We 
further  assume  a  nonnegative  IA  a(t)  =  1  +  Kb(t)  without  loss  of  generality.  The  narrowband 
constraint  can  then  be  summarized  similarly  as  in  [2]  via 


(26) 

K  «  1, 

(27) 

UJf  <C  CJc, 

(28) 

=  —  <  1, 

UJC 

(29) 

where  A  is  the  inverse  of  the  CR/FD,  which  is  also  called  the  FM  modulation  depth.  We  can 
easily  check  that  the  narrowband  constraint  requires  the  two  parameters  CR/FD  and  CR/IB  to  be 
sufficiently  large.  However,  it  is  not  the  case  in  the  wideband  scenario,  where  significant  errors 
are  incurred  regardless  of  the  choice  of  demodulation  approaches.  To  improve  the  demodulation 
performance  of  existing  approaches  in  the  wideband  scenario,  we  introduce  the  multirate  frequency 
transformations  (MFT)  approach  in  the  next  section,  that  increases  the  CR/FD  and  CR/IB  ratios 
for  FM  signals  and  for  AM-FM  signals  also. 


3  Methods,  Assumptions,  and  Procedures 

3. 1  Wideband  FM  Demodulation  Via  Multirate  F requency  Transformations 

The  MFT  approach  that  involves  multirate  systems  as  well  as  heterodyning  is  proposed  for  wide¬ 
band  FM  demodulation  that  usually  has  a:  1)  large  frequency  deviation;  2)  large  information  band¬ 
width.  The  proposed  approach  can  be  integrated  with  other  existing  FM  demodulation  methods, 
such  as  the  analytic  signal,  negative  feedback  demodulation  or  energy  operators,  serving  as  a  gen¬ 
eral  framework  for  wideband  FM  demodulation. 


3.1.1  Basic  MFT  System 

Using  the  scaling  property  of  the  Fourier  transform,  compression  in  frequency  domain  is  equivalent 
to  expansion  in  the  time  domain  expressed  as 

y(t)  =  x(at)  Y(u)  =  X  ,  (30) 

where  a  =  1  / R  <  1  is  the  factor  of  frequency  compression.  The  IF  of  the  compressed  signal 
becomes  a  scaled  version  of  the  input  IF  by  a  factor  R  expressed  as 


&i(t) 


u%(t) 

R 


UJr 


R+  RQi 


(31) 
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Note  that  for  a  compressed  signal,  the  carrier  frequency  is  also  scaled  by  the  same  factor  R,  which 
is  undesirable  since  the  ratios  CR/FD  and  CR/IB  that  we  wish  to  increase  still  remains  invariant. 
Hence  the  heterodyning  operator  is  cascaded  right  after  the  frequency  compression  module  in 
order  to  upshift  the  carrier  frequency  to  a  higher  level  where  we  can  attain  larger  CR/FD  and 
CR/IB  ratios.  The  compressed  signal  after  frequency  translation  and  bandpass  filtering  in  the 
heterodyning  module  is  given  by 


yUsh(t)  =  [ y  it)  cos  (udt)]  *  hBPF  (t) , 


(32) 


where  *  denotes  the  convolution,  ood  refers  to  the  amount  of  frequency  translation  and  hBpF (t) 
represents  the  impulse  response  of  the  bandpass  filter.  Specific  for  the  case  of  sinusoidal  FM,  it 
can  be  further  simplified  as 


Vushit)  ~  -A  cos 


(33) 


where  <f>(t )  denotes  the  phase  of  the  original  FM  signal.  The  resultant  signal  that  has  a  scaled 
information  bandwidth  with  a  higher  CR/IB  then  passes  through  the  demodulation  block  for  IF 
extraction.  Eventually,  the  IF  estimate  of  the  original  signal  is  evaluated  by  the  inverse  MFT 
relation  via 

u?*{t)  =  R(ui{m)-u}d),  (34) 


where  0Jl(t)  is  the  IF  of  the  compressed  and  frequency  translated  signal. 

As  for  a  discrete-time  signal,  compression  and  expansion  can  be  substituted  by  the  correspond¬ 
ing  multirate  operations  of  interpolation  and  decimation  respectively  as  described  in  [8]  with  their 
properties  carried  over  to  their  discrete  counterparts.  The  block  diagram  of  the  basic  MFT  frame¬ 
work  in  prior  work  [4]  is  illustrated  by  Fig.  1.  The  wideband  signal  is  first  sampled  above  the 
Nyquist  rate,  interpolated  by  a  factor  R  and  then  heterodyned  by  multiplication  with  cos(c odn), 
followed  by  a  discrete  Finite  Impulse  Response  (FIR)  bandpass  filter  with  a  scaling  module  based 
on  (33)  to  achieve  MFT.  Then  it  goes  through  a  demodulation  block  to  obtain  the  IF  estimation  of 
the  compressed  heterodyned  signal.  To  estimate  IF  of  the  original  signal,  the  compressed  hetero¬ 
dyned  IF  is  then  shifted  back  by  subtracting  ud,  decimated  by  R  and  scaled  back  approapriately 
according  to  (31),  followed  by  the  DAC  module. 


7 

Approved  for  public  release;  distribution  is  unlimited. 


Figure  1:  Block  diagram  of  the  basic  MFT  framework 


Interpolating  the  input  signal  will  result  in  the  reduction  of  both  the  frequency  deviation  and 
information  bandwidth  by  a  factor  of  R.  Similar  to  increasing  the  sampling  rate,  the  IF  of  the 
interpolated  signal  becomes  slow-varying  and  the  assumption  that  the  message  signal  remains 
constant  over  the  carrier  period  is  more  likely  to  hold,  which  in  turn  boost  the  performances  of 
conventional  demodulation  algorithms.  Meanwhile,  heterodyning  serves  the  purpose  of  increasing 
the  ratios  of  CR/FD  and  CR/IB  by  compensating  for  the  scaled  carrier  frequency.  By  passing 
through  the  heterodyning  process,  the  CR/FD  and  the  CR/IB  of  yush{t)  are  given  by 


CR' 

'CR 

Rud 

+ 

FD 

out 

[FD 

in  UJm 

CR 

'CR 

!  RuJd 

1B\ 

out 

IB  J 

in 

(35) 

(36) 


3.1.2  Error  Analysis  for  Wideband  FM  Signals 

Here  we  present  the  error  analysis  for  wideband  sinusoidal  FM  demodulation  in  ID  1 .  According 
to  the  handbook  of  mathematic  formula  and  basic  trigonometric  identities,  we  have 

OO 

cos(/3sin0)  =  J0(/3)  +  2^  J2fc(/3)  cos(2 kO), 

k= 1 

OO 

sin(/3  sin#)  =  2  E  J2k+1  ((3)  sin((2/c  +  1)6) , 

k= 0 

'The  partial  Hilbert  transform  is  a  separable  combination  of  the  ID  Hilbert  transforms,  consequently  this  analysis 
is  valid  for  the  proposed  2D  approach  as  well. 
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cos(w)  cos(v)  =  -[cos(w  —  v)  +  cos(u  +  n)], 


sin(w)  sin(u)  =  -[cos(w  —  v)  —  cos(u  +  v)], 

where  Jn(P)  denotes  the  Bessel  function  of  the  first  kind. 

Using  above  equations,  the  sinusoidal  FM  can  be  expanded  as 

cos(cuct  +  P  sin(u;mt)) 

=  cos(c oct)  cos  (P  sin(u;m£))  —  sin(u;cf)  sin  (P  sin(tcmf)) 


:  cos (uct)J0(P)  +  2  J2 k(P)  cos(2 kumt)  cos(o jct) 

k=L 

OO 

-  2  ^2  (P)  sin((2fc  +  1  )umt)  sin(tucf) 

k= 0 

oo 

:  cos (uct)J0(P)  +  2  ^2  'hk(P)  COS  ((uc  +  2 kum)t) 
k= 1 

OO 

+  2  J'2k (£)  cos  ((Wc  ~  2kuj  m)t) 

k=  1 
oo 

+  2  y]  J-2k+i(P)  cos  ((u)ct  +  (2k  +  l)t om)t) 

k= 0 
oo 

-  2  ^  J2fe+i(/0)  COS  ((wc  -  (2k  +  l)tcm)  t) 


k=0 


Its  quadrature  part  that  we  desire  for  demodulation  can  also  be  expanded  similarly  given  by 

sin  (ojct  +  P  sin(u;mf)) 

=  sin(tucf)  cos  (P  sin(cumf))  +  cos(cuct)  sin  (P  sm(ujmt)) 

OO 

=  sin (u>ct)J0(P)  +  2  J2jfc(/3)  cos(2 ku>mt)  sin(u;c£) 

k= 1 

OO 

+  2  ^2  J-2k+i  (P)  sin((2fc  +  1  )umt)  cos(cuct) 

k=0 

oo 

=  sm(u}ct)Jo(P)  +  2  ^2  J2 k(P)  sin  ((uc  +  2kum)t) 

k=  i 

OO  OO 

+  2  ^  J2fc(/5)  (sin(ccc  -  2kum)t)  +  2  ^  J2fe+i(/?)  sin  ((cucf  +  (2k  +  l)tcm)  f) 

/c=l  fc=0 

oo 

-  2  J2fe+i(/?)  sin  ((uc  -  (2k  +  l)o;m)  t) 

k= 0 
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Note  that  the  Hilbert  transform  of  the  corresponding  sinusoidal  FM  can  be  easily  obtained  by 
performing  the  Hilbert  transform  at  each  of  its  harmonics  as 

"H[cos  (u)ct  +  /3sin  (u)mt))] 

oo 

—  sm(cuct)Jo(/3)  +  2  ^  ^  ^2fc(/5)  sin  (c oc  H-  2/ccj  m)t ) 

k= 1 

oo 

+  2  ^2  J2 k{/3)  sin  (I uc  -  2kum\t ) 

fc=i 

OO 

+  2  ^  J2fc+i(/3)  sin  ((cuc  +  (2 k  +  l)ccm)  t) 

fc= o 

OO 

-  2^2  J2k+iC0)sin(|o;c  -  (2k  +  l)ujm\t) 

k= 0 

By  comparing  the  quadrature  and  the  Hilbert  transform  of  the  corresponding  sinusoidal  FM,  we 
observe  that  when  uc  <  num, 

"H[cos  ((cuc  —  nujm)t)]  =  sin (nuim  —  u>c)t, 

whereas  the  corresponding  part  in  the  quadrature  part  is  of  opposite  sign  sin(ccc  —  nujrn)t.  For  nar¬ 
rowband  sinusoidal  FM,  the  Bessel  coefficients  Jn(/3 )  associated  with  these  corresponding  side- 
lobes  are  very  small,  and  the  demodulation  error  incurred  by  the  Hilbert  transform  due  to  these 
opposite  sign  terms  can  be  neglected.  Hence  the  Hilbert  transform,  provides  almost  perfect  recon¬ 
struction  of  the  quadrature  part  we  desire  for  demodulation  using  the  analytic  signal,  for  sinusoidal 
FM,  in  the  narrowband  case.  For  extremely  wideband  sinusoidal  FM,  however,  due  to  its  large 
modulation  index  (3,  the  Bessel  coefficients  associated  with  these  sidelobes  cannot  be  neglected 
anymore.  As  a  result,  the  error  incurred  by  the  Hilbert  transform  is  much  more  significant.  In 
general,  this  unavoidable  error  imposes  a  lower  bound  on  the  error  performance  for  HT  and  energy 
operator  in  the  extremely  wideband  scenario. 

3.1.3  Alternative  MFT  System  for  Large  Conversion  Factors 

As  we  look  further  into  the  MFT  framework,  an  important  question  regarding  the  selection  of  the 
conversion  factor  R  arises.  Specifically  prior  work  [4]  only  deals  with  small  multirate  compression 
factors.  However,  larger  factors  over  hundred  or  thousand  can  be  supported  by  current  high-speed 
digital  signal  processing  with  large  memory.  It  is  intuitive  to  expect  a  further  reduction  in  the 
demodulation  error  since  the  gain  brought  by  frequency  compression  should  be  extendable  through 
the  use  of  a  larger  factor.  But  for  sufficiently  large  factors  R,  the  passband  of  the  lowpass  filter  in 
the  multirate  operation  and  that  of  the  heterodyne-band  pass  filter  (BPF)  operation  will  be  scaled 
by  R.  For  example,  if  R  =  1000,  we  require  a  lowpass  filter  with  cut-off  frequency  at  7r/1000 
and  a  bandpass  filter  with  a  passband  edge  less  than  or  equal  to  7r/1000.  However,  filters  with 
such  narrow  passbands  are  unrealistic  for  direct  implementation  by  any  structure2.  Therefore,  the 

2Narrow  passband  implies  clustered  poles  and  zeros  that  result  in  sensitivity  and  stability  issues  of  digital  filters  as 
described  in  [9], 
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design  of  the  BPF  within  the  basic  MFT  framework  becomes  the  bottleneck  that  limits  the  use  of 
a  very  large  factor. 

In  general,  the  MFT  framework  achieves  better  demodulation  performance  when  the  conver¬ 
sion  factor  becomes  larger.  However,  large  conversion  factors  require  the  heterodyning  bandpass 
filter  to  have  significantly  narrow  passband,  which  is  hardly  achieved  in  practice.  In  order  to  re¬ 
duce  the  burden  placed  on  the  practical  implementation  of  the  bandpass  filter,  we  first  consider  a 
different  MFT  framework  where  the  order  of  the  interpolation  operator  and  the  heterodyning  op¬ 
erator  are  exchanged.  Due  to  the  switch  of  interpolation  and  heterodyning,  the  CR/FD  and  CR/IB 
parameters  under  this  MFT  framework  are  given  by 


CR' 
FI)  _ 

out 

' CR ' 

FD 

ud 

+  — 

in 

(37) 

CR 

IB 

out 

'CR' 

IB 

U)d 

+  — . 

in  ^  f 

(38) 

By  comparing  these  ratios  with  (35)  and  (36),  note  that  the  upshift  frequency  cod  in  this  case  needs 
to  be  sufficiently  large  such  that  the  ratios  of  CR/FD  and  CR/IB  still  stay  at  high  level.  However, 
if  cud  is  too  large  such  that  the  resultant  carrier  frequency  after  heterodyning  exceeds  one  half  of 
the  sampling  rate,  then  we  have  to  interpolate  the  signal  first  by  an  appropriate  factor  in  order  to 
perform  discrete-time  bandpass  filtering  after  heterodyning.  Hence  the  practical  implementation 
of  MFT  framework  for  a  large  conversion  factor  is  not  as  simple  as  just  exchanging  the  order  of 
interpolation  and  heterodyning.  Actually,  an  interpolation  operation  is  still  required  prior  to  the 
heterodyning  with  an  appropriate  factor  that  depends  on  the  upshift  frequency  ud  and  the  sampling 
frequency  of  the  original  wideband  FM  signal.  This  implies  that  the  overall  interpolation  factor 
can  be  split  into  two  with  the  first  one  prior  to  the  frequency  translation  and  the  other  one  right 
after.  Then  upshifting  by  a  frequency  c od  that  is  not  too  large  would  result  in  a  relatively  small 
factor  for  the  first  interpolation,  thereby  lessening  the  burden  of  the  heterodyne-BPF. 

In  this  section,  we  propose  the  alternative  MFT  framework  for  large  conversion  factor  illus¬ 
trated  by  Fig.  2(a)  that  achieves  practical  design  of  the  heterodyne-BPF.  Where  Fig.  2(a)  is  the 
Block  diagram  of  the  proposed  MFT  framework  for  large  wideband  to  narrowband  conversion 
factors.  The  interpolation  module  of  the  prior  framework  is  separated  into  two  with  one  in  front 
of  the  heterodyning  module  and  the  other  one  right  after.  The  first  interpolation  module  has  a 
relatively  small  upsampling  rate  of  Ri  which  is  appropriately  chosen  such  that  the  discrete  BPF 
can  be  implemented  within  the  range  of  half  the  sampling  rate  after  heterodyning  the  signal  with 
a  frequency  translation  of  u>d.  The  relatively  small  R\  would  result  in  a  wider  passband  for  the 
discrete  bandpass  filter,  thus  reducing  the  its  design  of  complexity.  And  where  Fig.  2(b)  is  the 
Multistage  implementation  for  interpolation  modules  in  Fig.  2(a)  based  on  the  Noble  identity. 
Note  that  multistage  implementation  for  the  corresponding  decimation  modules  can  be  realized 
in  a  similar  fashion.  As  previously  discussed,  the  interpolation  module  of  the  prior  framework  is 
separated  into  two  with  one  in  front  of  the  heterodyning  module  and  the  other  one  right  after  in 
the  proposed  framework.  The  first  interpolation  module  has  a  relatively  small  upsampling  rate  of 
Ri  which  is  appropriately  chosen  such  that  the  discrete  BPF  can  be  implemented  within  the  range 
of  half  the  sampling  rate  after  heterodyning  the  signal  with  a  frequency  translation  of  cod.  The 
relatively  small  Ri  would  result  in  a  wider  passband  for  the  discrete  bandpass  filter,  thus  reducing 
complexity  of  the  design.  In  general,  there  is  a  sacrifice  in  terms  of  achievable  CR/IB  and  CR/FD 
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ratios  for  the  alternative  MFT  framework.  However,  the  basic  MFT  system  cannot  realize  large 
conversion  factors,  due  to  the  placement  of  impractical  constraints  on  the  BPF  design. 

In  addition,  a  binomial  smoothing  module  is  incorporated  into  the  alternative  framework  as 
shown  in  Fig.  2(a),  to  further  reduce  the  effects  of  noise.  Even  though  the  FM  signal  is  wideband, 
the  IF  waveform  itself  is  not  necessarily  wideband  in  nature.  In  many  cases,  the  wideband  FM 
is  primarily  generated  by  a  large  modulation  index  while  the  IF  waveform  still  remains  in  the 
narrowband  range.  Under  this  assumption,  by  applying  the  binomial  smoothing  we  can  efficiently 
filter  out  the  high  frequency  noise  in  the  corrupted  IF  estimation.  When  the  signal  to  noise  ratio 
(SNR)  is  high,  the  improvement  becomes  extremely  evident  as  we  shall  see  later.  Usually  we 
would  expect  a  gain  between  5  dB  and  10  dB  in  the  scenario  of  relatively  high  SNR. 


Heterodyning 


(a) 


(b) 


Figure  2:  Block  diagrams  of  the  alternative  MFT  system  for  large  conversion  factors  (a)  and  the 
Noble  Identity  relation  (b)  applied  in  the  system. 
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(a)  (b) 

Figure  3:  Comparison  between  performances  of  both  MFT  frameworks  under  noise  free  environ¬ 
ments.  (a)  Performance  of  the  basic  MFT  framework,  (b)  Performance  of  the  alternative  MFT 
framework. 


3.1.4  Examples  of  FM  Demodulation  via  Both  MFT  Frameworks 

In  this  section,  we  present  demodulation  results  using  the  basic  and  alternative  MFT  frameworks 
respectively  under  both  noise  free  and  noisy  environments.  Here  we  apply  the  GNGD  approach  as 
the  demodulation  technique  to  test  on  pure  FM  signals.  Note  that  the  demodulation  performance 
is  judged  by  the  normalized  RMS  IF  demodulation  error  (NRMSE)  throughout  this  section. 

We  first  look  at  the  case  of  a  wideband  sinusoidal  FM  signal  that  has  a  modulation  index  of 
10  and  the  CR/IB  of  20.  Under  a  noise  free  environment,  the  performance  of  the  basic  MFT 
framework  is  illustrated  by  Fig.  3(a).  Note  that  the  performance  associated  with  R  —  1,  i.e, 
the  origin  of  the  performance  curve  corresponds  to  direct  demodulation  by  GNGD  without  MFT, 
while  'jJd  is  the  normalized  upshift  radian  frequency  translation  in  the  range  of  [0, 7r].  By  applying 
a  large  conversion  factor  of  R  =  128,  a  reduction  of  around  40  dB  in  the  demodulation  error 
over  direct  GNGD  demodulation  is  attained.  The  result  of  Fig.  3(a)  confirms  the  claim  that  a 
large  conversion  factor  strengthens  the  benefits  achieved  by  frequency  compression  thus  leading  to 
significant  reduction  in  the  demodulation  error.  Moreover,  it  reflects  the  fact  that  the  use  of  a  larger 
factor  requires  a  very  high  order  FIR  bandpass  filter  with  a  satisfactory  frequency  response.  For 
example,  R  =  128,  demands  the  order  of  FIR  bandpass  filter  to  be  as  high  as  4096,  which  results  in 
unrealistic  parameters  for  the  narrow  passband.  This  constraint  seriously  limits  the  implementation 
in  a  practical  system  for  large  factors. 

To  relax  the  constraint,  the  alternative  MFT  framework  is  applied,  the  demodulation  perfor¬ 
mance  of  which  is  illustrated  by  3(b).  In  comparison  to  Fig.  3(a),  it  can  be  observed  that  the 
required  order  for  the  FIR  bandpass  filter  is  effectively  reduced  at  the  cost  of  sacrificing  a  small 
amount  of  improvement  in  demodulation  error.  For  instance,  when  R  =  128,  the  order  for  the  FIR 
bandpass  filter  drops  significantly  from  4096  to  512  with  just  4  dB  loss  in  performance,  suggesting 
no  observable  difference  in  performance  between  both  MFT  frameworks  except  for  the  dramatic 
reduction  of  the  order  for  the  heterodyne-BPF.  Note  that  the  frequency  response  of  the  heterodyne- 
BPF  in  the  alternative  MFT  framework  for  the  case  R  =  128  is  illustrated  by  Fig.  4(b),  which  has 
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Figure  4:  Comparison  between  the  basic  and  alternative  MFT  frameworks  with  a  conversion  factor 
R  =  128.  (a)  Previous  MFT  framework,  (b)  Alternate  MFT  framework. 


much  wider  passband  compared  to  the  frequency  response  of  the  BPF  in  the  basic  MFT  framework 
shown  in  Fig.  4(a)  and  thus  practical  for  implementation. 

For  noisy  environments,  the  performance  of  the  basic  MFT  framework  is  illustrated  by  Fig.  5(a) 
for  the  same  sinusoidal  FM  signal  corrupted  by  additive  white  Gaussian  noise  (AWGN).  Based  on 
this  observation,  improvement  over  GNGD  demodulation  alone  varies  according  to  the  different 
SNRs.  For  example,  when  the  SNR  is  20  dB,  the  improvement  is  only  aoumd  6  dB,  and  when 
the  SNR  increases  to  40  dB,  improvement  increases  to  around  20  dB.  In  addition,  Fig.  5(b)  sum¬ 
marizes  the  results  of  the  alternative  MFT  framework  in  the  presence  of  noise,  indicating  better 
performance  due  to  the  binomial  smoothing.  For  the  case  of  20  dB  SNR,  the  improvement  in¬ 
creases  to  15  dB  compared  with  the  basic  MFT  system.  Note  that  the  NRMSE  gradually  becomes 
saturated  as  R  increases,  due  to  LTI  filtering  induced  harmonic  distortion  of  the  FM  signals3.  In 
contrast,  the  proposed  MFT  approach  results  in  a  better  performance,  and  offers  a  more  loose 
constraint  on  the  BPF  design. 

We  further  investigate  an  extreme  wideband  scenario  under  the  noise  free  environment,  where 
the  modulation  index  /3  is  as  large  as  50  and  the  frequency  deviation  is  equal  to  the  carrier  fre¬ 
quency  with  the  IF  varying  over  the  entire  carrier  range.  For  the  signal  of  interest,  the  IF  estimates 
of  both  the  basic  and  alternative  MFT  frameworks  are  illustrated  by  Fig.  6(a)  and  Fig.  6(b)  re¬ 
spectively.  It  can  be  observed  that  the  GNGD  demodulation  alone  fails  in  this  extreme  wideband 
scenario,  while  both  MFT  frameworks  maintain  tracking.  The  observation  implies  that  both  MFT 
frameworks  guarantee  the  demodulation  with  acceptable  performance  even  in  the  worst  scenarios 
where  conventional  algorithms  would  normally  fail. 

To  quantify  the  performance  of  the  MFT  approach,  we  explore  another  scenario  where  the  sig¬ 
nal  is  a  wideband  linear  chirp  instead  of  a  sinusoidal  FM.  The  short-time  spectrum  of  the  chirp  sig¬ 
nal  is  illustrated  by  Fig.  7(a).  To  validate  the  performance  of  the  MFT  approach  ,  we  can  compare 

'Since  they  are  only  approximate  eigenfunctions  as  described  in  [10]. 
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Figure  5:  Comparison  between  performances  of  both  MFT  frameworks  in  environments  corrupted 
by  AWGN.  (a)  The  basic  MFT  framework,  (b)  The  alternative  MFT  framework. 


*  P=50,CR/IB=50,R=1 28,L=4096,w  =0.31 41 6 

rr  K  d 


(a) 


(b) 


Figure  6:  Demodulation  performances  of  both  MFT  frameworks  with  conversion  factor  R  =  128 
and  normalized  radian  frequency  shift  w<i  =  0.17T  under  the  extreme  senario  with  modulation  index 

P  =  50. 


the  variance  of  error  with  respect  to  the  chirp  rate  estimate  with  its  Cram  'er  Rao  lower  bound 
(CRLB). 

The  chirp  rate  can  be  obtained  from  the  demodulated  IF  followed  by  a  least  square  estimator.  In 
the  presence  of  noise  with  different  SNRs,  the  result  is  summarized  in  Fig.  7(b).  Improvement 
over  the  GNGD  alone  is  more  apparent  with  respect  to  lower  SNR,  indicating  satisfactory 
performance  of  the  MFT  approach  in  noisy  conditions.  Also  note  that  gap  between  the  error 
variance  estimate  of  the  MFT  approach  and  the  corresponding  CRLB  is  nearly  a  constant,  which 
can  be  explained  directly  via  the  loss  of  spectrum  incurred  due  to  filtering  of  the  FM  signal. 
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Short-time  Spectrum  of  the  chirp  signal 


Time  Index  SNR(dB) 

(a)  (b) 


Figure  7:  Performance  of  the  alternative  MFT  framework  with  a  multirate  conversion  factor  of  128 
in  wideband  linear  chirp  scenario. 


3.1.5  MFT  Compatibility  With  HT  And  ESA  Approaches 

In  general,  the  MFT  framework  can  be  combined  with  most  demodulation  approaches  to  enhance 
their  performances  in  the  wideband  scenario.  But  the  gain  achieved  by  using  the  MFT  framework 
varies  with  the  demodulation  algorithm  employed.  For  a  given  demodulation  scheme,  the  gain  in 
error  reduction  via  the  Bi-directional  multirate  frequency  transformation  (BMFT)  approach  mainly 
depends  on  its  inherent  sensitivity  to  CR/IB,  CR/FD,  and  the  modulation  index  of  the  image.  Here 
we  refer  to  the  comparison  results  of  the  HT  approaches  and  ESA  approaches  in  ID  [1], 

The  HT  demodulation  approaches  are  usually  not  very  sensitive  to  the  CR/IB  of  the  signal.  In 
general,  the  demodulation  performance  of  the  Hilbert  transform  is  invariant  to  the  CR/IB  of  the 
signal.  As  shown  in  the  work  [1]  by  Potamianos  and  Maragos,  the  error  vesus  CR/IB  response 
of  the  Hilbert  transform  is  almost  flat,  suggesting  that  by  increasing  the  CR/IB  via  the  MFT,  the 
gain  we  achieve  through  error  reduction  is  very  limited.  On  the  other  hand,  the  ESA  is  shown  to 
be  very  sensitive  to  the  CR/IB  of  the  signal,  where  the  error  decreases  sharply  as  the  CR/IB  of  the 
signal  increases.  Hence  the  gain  attained  by  increasing  the  CR/IB  via  the  MFT  with  demodulation 
schemes  based  on  ESA  is  expected  to  be  more  significant  than  demodulation  schemes  using  the 
HT  approaches. 

As  for  the  CR/FD,  both  the  HT  and  ESA  demodulation  methods  suffer  from  large  deviation 
compared  with  the  carrier  frequency.  Hence  the  BMFT  becomes  beneficial  to  both  of  them  via  an 
increase  in  the  CR/FD  of  the  input. 

Since  the  modulation  index  is  jointly  determined  by  both  the  CR/IB  and  CR/FD,  its  influence 
on  the  demodulation  performance  is  not  straightforward.  However,  the  modulation  index  of  the  FM 
signal  cannot  be  too  large  such  that  the  spectrum  of  the  sidelobes  around  the  origin  are  significant 
and  begin  to  incur  significant  demodulation  error.  The  error  incurred  due  to  a  large  modulation  in¬ 
dex  is  unavoidable  for  any  demodulation  algorithm.  It  imposes  a  lower  bound  on  the  demodulation 
error.  Even  the  MFT  framework  cannot  reduce  the  error  significantly  in  this  extreme  wideband 
scenario. 
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In  conclusion,  the  proposed  MFT  approach  is  more  suitable  for  the  ESA  demodulation  schemes. 
Moreover,  the  MFT  framework  helps  to  overcome  the  estimation  range  constraint  inherent  in  the 
ESA  demodulation  schemes.  By  compressing  and  shifting  the  original  input  in  the  frequency  do¬ 
main,  we  can  always  choose  appropriate  conversion  factors  and  translation  frequencies  to  convert 
its  frequency  components  into  the  range  between  0  and  |.  This  is  another  advantage  for  combining 
the  MFT  framework  with  the  ESA  demodulation  schemes. 

In  addition,  the  MFT  framework  can  be  extended  in  a  straightforward  separable  way  into  multi¬ 
dimensional  space  and  combined  with  the  multi-dimensional  higher  order  differential  operators 
recently  proposed  by  Salzenstein  et  al.  [11].  Hence  the  analysis  above  also  applies  to  the  BMFT 
framework  introduced  in  the  next  chapter. 
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4 


Results  and  Discussions 


4.1  Wideband  Image  Demodulation 

The  amplitude-modulation  frequency-modulation  (AM-FM)  model  has  found  various  applications 
with  images  recently  including  image  analysis,  texture  processing  and  fingerprint  classification 
[12].  According  to  earlier  work  [13]  by  Havlicek,  Bovik  et  al.,  nonstationary  images  can  be  mod¬ 
eled  as  superpositions  of  multiple  AM-FM  components. 


n 

I(x,y)  =  ^2  ai(x,y)  cos  (<j>i(x,y)) .  (39) 

2=1 


The  multi-component  AM-FM  image  is  first  decomposed  by  employing  a  set  of  bandpass  filters 
such  as  Gabor  filterbanks  or  via  the  use  of  the  Bi-dimensional  Empirical  Mode  Decomposition 
(BEMD).  Each  resulting  monocomponent  AM-FM  image  is  further  demodulated  into  correspond¬ 
ing  instantaneous  amplitude  (LA)  a(x,  y )  and  instantaneous  frequency  vector  (IF), 


V<t>(x,y) 


d <f>(x,y)  d <t>{x,y) 


T 


dx  ’  dy 


(40) 


In  particular,  the  IA  depicts  the  contrast  present  in  the  image,  while  the  IF  reveals  the  locally 
emergent  frequency  variation.  Conventional  image  demodulation  approaches  involve  2D  extension 
of  the  analytic  signal  (AS)  and  multidimensional  Energy  Separation  Algorithm  (ESA)  [14]  with 
additional  processing  techniques  such  as  Dominant  Component  Analysis  (DCA).  However,  in  most 
of  these  approaches,  narrowband  assumptions  were  imposed  on  each  AM-FM  component  of  the 
image.  For  example,  most  literature  implicitly  assumes  the  AM-FM  image  to  be  globally  wideband 
yet  each  of  its  components  to  be  locally  narrowband.  In  general,  both  the  IA,  a(x,  y)  and  the  IF, 
V<j)(x,y)  of  a  single  component  are  assumed  to  be  slowly  varying,  otherwise  the  approximations 
inherent  in  most  demodulation  approaches  are  no  longer  valid  and  incur  significant  error  especially 
under  the  wideband  scenario. 

In  this  chapter,  we  propose  the  2D  extension  of  the  MFT  approach  called  the  Bi-dimensional 
Multirate  Frequency  Transformations  (BMFT)  that  can  be  combined  with  a  variety  of  demodula¬ 
tion  techniques  to  enhance  the  demodulation  performance,  traditionally  limited  by  the  narrowband 
constraint  on  the  frequency  modulation  part  of  the  monocomponent  AM-FM  image. 


4.1.1  Wideband  Frequency  Modulated  Image 

Since  images  are  not  separable  in  general,  we  cannot  directly  extend  the  related  definitions  of  ID 
wideband  signal  to  2D.  For  instance,  we  are  not  able  to  define  the  corresponding  modulation  index 
along  any  specific  direction  for  2D  images  in  a  global  sense.  However,  as  proposed  by  Pattichis 
and  Bovik  [15],  the  complex  FM  image  can  be  locally  approximated  by  the  product  of  two  ID 
FM  signals.  The  corresponding  ID  signals  are  defined  along  the  directions  of  the  eigenvectors 
of  the  instantaneous  frequency  gradient  tensor  (IFGT),  which  is  simply  the  Hessian  of  the  phase. 
Fet  z  —  [z] ,  z2\  denote  the  representation  under  the  eigenvector  coordinate  system,  <p(z)  denote 

the  phase  of  the  FM  image  and  F  denote  the  IFGT.  According  to  [15],  around  a  given  point 
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z0  =  [ai,  a2\  of  the  image,  the  local  phase  is  approximated  using  a  Taylor  series  expansion  by 

4>{z)  ~  4>{a i,  a2)  +  (f>i(zi)  +  02(z2),  (41) 


where: 


7/  \  d<t> ,  w  x  .  Aj(ai,a2), 

=  ^7(01,  a-2){zi  -  at)  H - - - (2j 


aj 


i  =  1,2. 


(42) 


Note  that  Ai  and  A 9  are  the  eigenvalues  of  F.  The  complex  FM  image  is  then  locally  approxi¬ 
mated  by  the  product  of  two  ID  FM  functions  defined  with  respect  to  the  eigenvector  coordinate 
directions: 


exp  \j</>(zi,  z2)\  «  exp  \j<pi{zi)}  exp  \jfa(z2)]. 


(43) 


As  a  result,  FM  images  can  be  viewed  as  locally  separable,  and  we  can  define  the  wideband  FM 
images  locally  by  deriving  the  corresponding  parameters  along  the  eigenvector  coordinate  direc¬ 
tions.  For  example,  if  the  modulation  index  along  the  directions  of  the  IFGT  eigenvectors  are 
sufficiently  large  around  the  reference  point,  the  image  can  be  expected  to  exhibit  similar  wide¬ 
band  pattern  locally  analogous  to  the  ID  case.  Assume  the  modulation  index  along  both  directions 
of  the  eigenvectors  are  denoted  by  f3Zl  and  JZ2  respectively.  To  locally  indicate  whether  an  FM 
image  is  wideband  or  not  in  the  sense  similar  to  the  literature  in  FM  communication  systems,  we 
can  define  a  local  modulation  index  (3i  via 


A  =  (/*!,+/£) 1/2  • 


(44) 


If  (5i  is  sufficiently  large,  then  the  corresponding  FM  image  exhibits  a  wideband  pattern  locally. 
This  definition  of  the  local  modulation  index  (3i  is  especially  effective  for  the  sinusoidal  FM  images 
with  IF  components  specified  along  the  horizontal  and  vertical  directions.  In  this  case,  the  modula¬ 
tion  index  for  each  eigenvector  direction  can  be  well  defined  according  to  the  original  formulation 
in  the  ID  context. 


4.1.2  Partial  Hilbert  Transform  Demodulation 

The  partial  or  directional  Hilbert  transform  based  on  analytic  image  is  widely  used  for  monocom¬ 
ponent  AM-FM  image  demodulation. 

It  can  be  derived  simply  from  the  ID  Hilbert  transform.  Let  f(x)  :  Mn  — »  M  and  et  denote 
the  unit  vector  in  the  Xi  direction.  The  partial  Hilbert  transform  along  the  direction  of  e,  is  then 
defined  by 

[/W]  =  [  (45) 

Jr 

Note  that  the  partial  Hilbert  transform  defined  here  is  specified  along  the  horizontal  and  verti¬ 
cal  directions.  Thus  it  is  particularly  suitable  for  demodulating  images  with  frequency  modulation 
patterns  emerging  along  the  same  directions.  Assume  that  I(x,y )  is  a  monocomponent  AM-FM 
image  with  frequency  modulation  patterns  emerging  along  both  the  horizontal  and  vertical  direc¬ 
tions  given  by 

I(x,  y )  =  A{x,  y)  cos  +  Slyy  +  (p(x,  y)) ,  (46) 
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where  f lx  and  ily  denote  the  carrier  (mean)  frequencies  along  the  x  and  y  directions  respectively. 
Under  this  assumption,  analogous  to  the  one  dimensional  problem,  we  are  able  to  approximate  the 
partial  Hilbert  transform  along  the  x  axis  by  fixing  variable  y  via 

Hx  [I (x,  y)}  «  A(x,  y)  sin  (Qxx  +  (Qyy  +  ip(x,  y))) .  (47) 


This  approximation  holds  only  under  certain  conditions  similar  to  the  ID  case.  First  of  all,  the 
A(x,y)  should  be  slow-varying  and  narrowband.  Analogous  to  the  one  dimensional  narrowband 
constraint,  Qx  is  required  to  be  sufficiently  large  such  that 


> 


d<p(x,y) 

dx 


(48) 


A  similar  condition  applies  to  the  y  axis  as  well.  As  a  result,  we  can  easily  obtain  the  approx¬ 
imations  for  the  IA  and  the  IF  components  with  respect  to  both  directions  according  to  the  2D 
extensions  of  the  corresponding  one  dimensional  IA  and  IF  estimates. 

Synthetic  examples  are  provided  to  illustrate  the  efficacy  of  the  proposed  approach.  Since  the 
synthetic  images  we  present  are  predefined  with  its  frequency  modulation  patterns  along  horizon¬ 
tal  and  vertical  directions,  the  partial  Hilbert  transform  is  preferred  over  the  monogenic  image 
approach  to  serve  as  the  comparison  technique  to  the  proposed  approach  due  to  the  following  rea¬ 
sons:  1)  the  demodulation  performance  of  the  partial  Hilbert  transform  is  as  competitive  as  the 
Riez  transform  when  the  image  display  frequency  modulation  patterns  along  the  same  specified 
directions,  2)  the  partial  Hilbert  transform  can  be  implemented  easily  via  the  ID  Hilbert  transform. 


4.1.3  Higher  Order  Energy  Operator 


A  variety  of  methods  based  on  the  multidimensional  energy  operator  [14]  are  also  widely  used  for 
AM-FM  image  demodulation.  An  image  demodulation  algorithm  based  on  higher  order  Teager- 
Kaiser  operators  is  proposed  in  the  recent  work  [16]  by  Salzenstein,  Diop  and  Boudraa.  They  have 
been  reported  to  provide  better  performance  for  narrowband  AM-FM  images  than  the  classical  2D 
ESA  [14]. 


by 


The  k- order  differential  energy  operator  (DEO)  [20]  in  ID  for  any  given  signal  s(t)  is  defined 


Vk[s(t)] 


ds(t)  dk  1s(t)  ,  d ks(t) 

dt  dtk~l  dtk 


(49) 


where  \P2  refers  to  the  commonly  used  Teager-Kaiser  energy  operator.  For  a  given  image  I (k.  /), 
the  discrete-time  higher  order  demodulation  algorithm  (DHODA)  [16]  can  be  summarized  via: 
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I1(k,l)  =  1-[I(k+l,l)-I(k-l,l)], 
I2(k,l)=1-[I(k,l+l)-I(k,l-l)], 
h2(k,l)=l-{I2(k+l,l)-I2(k-l,l)], 
f2[-f  (k,  ()]  ={2  [I(k,  i)]2  -  I(k  -  1  ,l)I(k  +  l,l)-I(k,l-  1  )I(k,  l  +  1)} 

+  2[/,(fc,  l)I2(k,l)  -  I(k,  l)In(k,l)], 

I}2(k,l)  =  llh2(k  +  1,1)  -  I12(k  -  1,1)], 
Il2(k,l)=1-[h2(k,l  +  l)-h2(k,l-l)\, 


/*2[A(fc,0]»2[  I2(k,l)]\1/2 
\H2{I}2(k,l)  +  lUk,l)]J 


|«i(*,oi 


=  arcsin 


/WJ 

WW)] 


(50) 

(51) 

(52) 

(53) 

(54) 

(55) 

(56) 

(57) 


fe(M)  I 


=  arcsin 


('Hh2(k,l)] 

v  «'2[n  (M)] 


(58) 


where  a(k,  1)  is  the  IA  estimation  while  Vli(k,  l )  and  Q^fA-,  /)  are  the  IF  estimations  along  the  spa¬ 
tial  axes  of  the  image.  Note  that  the  IF  estimations  are  obtained  through  the  inverse  sine  function, 
indicating  that  the  demodulation  approaches  based  on  the  ESA  can  only  estimate  IF  components 
that  range  between  0  and  or  in  other  words,  up  to  one  fourth  of  the  sampling  frequency.  More¬ 
over,  the  demodulation  approaches  based  on  the  ESA  also  suffer  from  the  narrowband  constraint 
as  the  Hilbert  transform.  Both  the  IA  and  the  IF  waveforms  may  not  vary  too  fast  or  too  greatly  in 
value. 


4.1.4  Bi-dimensional  Multirate  Frequency  Transformations 

The  Multirate  Frequency  Transformations  (MFT)  approach  that  performs  wideband  to  narrowband 
conversion  of  a  FM  signal  was  proposed  in  the  previous  section.  The  MFT  primarily  increases  the 
CR/IB  and  CR/FD  ratios  of  the  original  signal  to  improve  the  demodulation  performance  of  the 
conventional  demodulation  approaches.  The  corresponding  Bi-dimensional  Multirate  Frequency 
Transformation  (BMFT)  is  derived  by  generalizing  the  underlying  idea  to  2D.  For  simplicity,  as¬ 
sume  that  the  input  is  a  monocomponent  wideband  FM  image  of  the  form 

J(x,y)  =  A  cos  (</>(x,y)) .  (59) 

It  is  first  compressed  in  frequency  domain  by  appropriate  factors  R  =  diag|/Tc,  Ry\,  which  corre¬ 
sponds  to  spatial  expansion  given  by 

Ji(x,y)  =  2lcos  ^0  •  (60) 
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Then  we  heterodyne  the  resultant  image  by  a  frequency  translation  vector  D  =  Qr,  Qy]  via 


J2(x,  y )  =Ji(x,  y )  cos(Oxx)  cos(fiyy) 


A 


x  y 


=  — cos  0^  +  0  — cos(f lyy) 


Rt  R, 


A 

+  -COS 


A 


ttxx -0  (  TT,^-)  )  COS (Qyy) 
\  ti'X  /  / 


=  - cos  nxx  +  Qyy  +  0  (  -5-,  -5- 


x  y 


Rx  Ry 


A 


+  -  cos  nxx  -  nyy  +  0  (  -5-,  -5- 


x  y 


A 


Rx  Ry 
x  y 


+  —  cos  +  Vtyy  -  0  ,  — 

+  j  cos  ^Qxx  -nyy-<f>  (j^,  JJ- 


(61) 


After  the  heterodyning  a  2D  bandpass  filter  is  applied  to  extract  the  desired  high-frequency  term 
through 


J (x,  y)  =  J2(x,  y)  *  hBP(x,  y) 


A 

~  —  cos 
4 

A 

=  —  cos 
4 


^flxx  +  Vtyy  +  0 
(0(24  2/))  • 


(62) 


Assume  that  the  support  of  the  original  image  spectrum  J(u\,  u2)  is  within  the  range  ay  e  |— Qt,  Q,], 
1,2.  The  2D  bandpass  filter  should  be  carefully  designed  with  its  passband  range  given  by 


|tCi|  £  [D-e,  O-e  T 


Vll. 

R, 


|^2  |  £  [^j/,  Qy  + 
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(63) 


On  one  hand,  the  carrier  (or  mean)  frequencies  of  the  modulation  in  both  dimensions  are  increased 
via  the  frequency  translation  vector  Q  =  [Ox,  Q0.  On  the  other  hand,  the  bandwidth  of  the  mod¬ 
ulating  image  is  reduced  by  the  appropriate  conversion  factors  R  =  diag [Rx,Ry].  These  two 
benefits  derived  via  the  BMFT  are  crucial  to  improving  IF  demodulation  due  to  the  following 
reasons:  1)  a  majority  of  demodulation  approaches  require  the  input  to  have  high  CR/IB  in  both 
dimensions,  2)  the  CR/FD  of  the  input  must  be  sufficiently  large  in  both  dimensions  such  that  the 
deviations  of  the  IF  components  can  be  tolerated  in  accordance  with  the  narrowband  assumption 
of  the  input. 

Then  we  can  recover  the  IF  of  the  input  image  from  the  IF  estimation  of  the  transformed  image 
J(x,  y)  via  the  inverse  BMFT.  Assume  that  the  IF  components  of  -J(X,  y)  and  J(x,  y)  are  given  by 


d(p(x,y) 

dx 


^2  (x,y) 


d(p(x,y) 

dy 


(64) 
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(65) 


«■(*•»)  =  ^-  «»(*•»)  =  ^- 

The  IF  estimation  for  input  image  is  recovered  through  the  inverse  BMFT  relation  (see  Fig.  8) 


via 


^1  (•£?  2/)  Fix  (  (FixX)  Ry'y')  Vla 


(66) 


D2(z,  y )  =  Ry  (n2(Rxx,  Ryy )  -  ,  (67) 

where  Cli(Rxx,  Ryy )  and  Q2(RxX,  Ryy)  represent  spatial  compression  (or  frequency  expansion) 
of  the  IF  estimation  for  the  transformed  image  J(x,y). 

In  order  to  implement  the  BMFT  in  discrete  time,  we  replace  compression  and  expansion  in 
frequency  domain  by  their  discrete  equivalences.  Note  that  the  compression  in  the  frequency  do¬ 
main  corresponds  to  interpolation  while  the  expansion  corresponds  to  decimation.  As  a  result, 
the  block  diagram  of  the  BMFT  demodulation  approach  is  depicted  in  Fig  8.  The  BMFT  is  im¬ 
plemented  through  discrete-time  operations  of  interpolation,  heterodyning,  and  bandpass  filtering. 
The  transformed  FM  image  is  then  demodulated  via  a  monocomponent  demodulation  approach 
and  the  original  IF  components  are  recovered  from  the  IF  estimation  of  the  transformed  image  via 
the  inverse  BMFT.  However,  discrete-time  operations  in  2D  may  not  be  extended  in  a  straightfor¬ 
ward  way  from  their  ID  counterpart.  It  primarily  depends  on  whether  the  input  image  is  separable 
or  not.  If  the  image  is  separable,  each  operation  of  the  BMFT  can  be  implemented  by  simply 
cascading  its  ID  equivalence  with  action  along  each  dimension.  If  this  is  not  the  case,  we  need 
to  pay  attention  to  a  few  issues  associated  with  the  non-separable  nature  of  the  image,  which  are 
often  intractable.  Therefore,  deriving  2D  operations  directly  from  their  ID  realizations  is  favored 
in  terms  of  its  simplicity  for  practical  implementation. 

In  particular,  the  frequency  upshift  in  heterodyning  is  not  uniquely  defined  in  2D.  This  can 
be  realized  either  through  the  product  of  two  separable  cosine  terms  cos  (Hi m)  cos  (02n)  in  each 
dimension  as  in  Fig.  8  or  via  just  one  cosine  term  cos  (Dim  +  02n)  in  the  diagonal  direction. 
These  two  different  realizations  will  lead  to  different  designs  of  the  2D  bandpass  filters.  For  the  first 
case,  we  only  need  to  bandpass  the  outermost  quarter  of  the  frequency  spectrum  in  each  quadrant 
to  extract  the  desired  high  frequency  term  expressed  in  (62).  This  is  simply  achieved  using  a 
separable  2D  bandpass  filter,  by  cascading  two  ID  bandpass  filters  with  action  in  each  dimension. 
For  the  latter  case,  where  the  frequency  upshift  is  diagonal,  it  is  much  more  complicated  to  achieve 
the  same  objective  with  a  realizable  design  of  the  2D  bandpass  filter.  The  issue  is  that  the  2D 
frequency  upshift  in  the  diagonal  direction  results  in  only  two  copies  of  the  original  spectrum 
in  the  first  and  third  quadrant  respectively.  However,  the  design  of  a  2D  bandpass  filter  whose 
passband  is  only  present  in  two  quadrants  and  also  capable  of  extracting  the  desired  term  in  (62)  is 
computationally  complex. 

In  addition,  the  BMFT  framework  can  be  directly  applied  on  a  monocomponent  AM-FM  image 
provided  that  the  IA  of  the  given  image  is  slowly  varying,  which  is  inherently  assumed  by  most 
image  demodulation  algorithms.  Under  this  constraint,  the  IA  can  be  approximated  via  the  inverse 
BMFT  relation  by 

A(x,  y)  =  A  (Rxx,  Ryy ) ,  (68) 

where  A  (. x ,  y)  denotes  the  demodulated  IA  of  the  transformed  image  J(x.  y).  Numerical  results 
are  given  later  in  this  chapter  to  support  this  claim. 
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Figure  8:  Block  diagram  of  the  Bi-dimensional  Multirate  Frequency  Transformations  (BMFT). 


4.1.5  Example  of  Wideband  Image  Demodulation 


In  this  section,  we  present  numerical  results  for  the  proposed  BMFT-energy  approach.  We  begin 
with  an  example  of  a  synthetic  sinusoidal  AM-FM  image  as  illustrated  in  Fig.  9.  The  expression 
of  this  synthetic  image  and  its  corresponding  IF  components  are  given  by 


f(m,  n)  = 


(  7T  7T  \1  7T  77  f  7T  77\ 

1  +  0.5  cos  (  — m  H - n  cos(  —m  H — n  +  6  sin  — m  H —  + 
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V50 
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/  7 T 
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77  \\ 
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„  ,  .  d<f)(m,n)  77  77  (  77  \ 

VL2(m,  n)  =  — - =  — I —  cos  (  — n )  .  (71) 

v  ’  ;  dn  3  6  V30  / 

Note  that  the  frequency  modulation  index  is  6  along  the  horizontal  direction  and  5  along  the 
vertical  direction,  both  of  which  are  sufficiently  large.  We  can  also  easily  check  that  the  CR/IB  and 
the  CR/FD  are  small  along  both  directions  for  this  wideband  sinusoidal  example. 

The  demodulation  results  via  the  DHODA  combined  with  the  proposed  BMFT  framework 
using  conversion  factors  of  R  —  diag[8,  8]  are  given  by  Fig.  10.  Specifically,  Fig.  10  represents  the 
perspective  plot  of  the  demodulation  results  for  the  wideband  sinusoidal  example  via  the  DHODA 
and  via  the  BMFT-DHODA  with  conversion  factors  R  =  diag[8,  8].  Fig.  10(a)  is  the  estimation  of 
the  IA  via  the  DHODA.  Fig.  10(b)  is  the  estimation  of  the  IA  via  the  BMFT-DHODA.  Fig.  10(c) 
is  the  estimation  of  the  IF  component  along  the  horizontal  direction  via  the  DHODA.  Fig.  10(d) 
is  the  estimation  of  the  IF  component  along  the  vertical  direction  via  the  DHODA.  Fig.  10(e)  is 
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Figure  9:  Synthetic  sinusoidal  AM-FM  image  with  wideband  FM  components. 


the  estimation  of  the  IF  component  along  the  horizontal  direction  via  the  BMFT-DHODA.  Finally, 
Fig.  (f)  10(f)  is  the  estimation  of  the  IF  component  along  the  vertical  direction  via  the  BMFT- 
DHODA.  It  is  compared  with  the  demodulation  via  the  DHODA  alone.  Note  that  the  demodulated 
IF  with  respect  to  either  the  horizontal  or  the  vertical  direction  exhibits  a  sinusoidal  pattern  along 
that  direction,  which  can  be  inferred  from  (69).  We  can  easily  observe  that  both  the  IA  and  the  IF 
obtained  via  the  BMFT-DHODA  are  smoother  than  that  via  the  DHODA  alone.  The  MSE  for  the 
IA  is  reduced  from  0.8910  to  0.0649  and  the  root  mean  square  error  (RMSE)  for  the  IF  is  reduced 
significantly  from  1.83%  to  0.25%  through  the  BMFT  framework.  Here  we  define  the  RMSE  as 
the  1,2  norm  of  the  difference  between  the  true  IF,  V(f>(x,  y )  and  the  estimated  IF,  Vo{x.  y )  against 
the  l2  norm  of  the  true  IF  itself  via 


RMSE 


II  V0(%,  y)  -  vfer/)  Ik 

l|V0(x,j/)||l2 


x  100%. 


(72) 


In  fact,  the  demodulation  error  can  be  further  reduced  if  larger  multirate  conversion  factors 
are  applied,  similar  to  the  ID  case  of  MFT.  As  shown  in  Table  1,  the  demodulation  errors  of  the 
proposed  BMFT  framework  using  DHODA  with  conversion  factors  of  R  —  diag[8,  8]  (BMFT- 
DHODA-8)  and  with  conversion  factors  of  R  —  diag[16, 16]  (BMFT-DHODA- 16)  are  compared 
with  the  DHODA  alone.  The  RMSE  can  be  further  reduced  to  0.18%  via  the  BMFT  with  con¬ 
version  factors  R  =  diag[16.16],  achieving  an  error  reduction  of  10  times  compared  with  the 
DHODA  alone.  The  use  of  larger  factors  results  in  a  much  narrower  passband  for  digital  filters  that 
are  difficult  to  realize  in  practice  due  to  a  sharper  transition  band.  To  implement  such  FIR  filters, 
a  very  large  filter  order  is  required,  which  may  not  be  acceptable  in  terms  of  the  desired  system 
complexity.  Therefore  the  choice  of  multirate  conversion  factors  should  be  weighted  between  the 
tolerance  for  demodulation  error  and  that  of  the  system  complexity.  Recently  the  BMFT  structure 
using  large  conversion  factors  has  been  implemented  via  an  equivalent  structure  that  alleviates  the 
complexity  of  practical  realization.  Aided  by  this  equivalent  structure,  we  may  further  reduce  the 
demodulation  error  by  exploring  even  larger  conversion  factors. 

The  demodulation  results  via  the  partial  Hilbert  transform  alone  (HTDA)  and  the  BMFT  and 
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(C) 


(d) 


(e)  (f) 

Figure  10:  Perspective  plot  of  the  demodulation  results  for  the  wideband  sinusoidal  example  via 
the  DHODA  and  via  the  BMFT-DHODA 


the  partial  Hilbert  transform  combination  with  conversion  factors  R  =  diag[8,  8]  (BMFT-HTDA-8) 
are  also  illustrated  in  Table  1.  We  can  observe  that  the  error  reduction  via  the  BMFT  framework  for 
the  HT  demodulation  is  not  as  obvious  as  for  the  ESA  demodulation,  which  justifies  our  discussion 
in  the  previous  section. 

In  the  second  experiment,  we  test  the  BMFT  and  the  DHODA  combination  on  a  wideband 
sinusoidal  FM  image  whose  IF  components  are  out  of  the  range  [0,  |],  which  is  the  estimation 
range  constraint  for  demodulation  approaches  based  on  the  ESA.  As  shown  in  Fig.  11,  the  synthetic 
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Table  1 :  Comparison  of  the  demodulation  errors 


var(^)  MSE(a,a)  var(^=^)  var(^k) 


DHODA 

4.7876 

0.8910 

2.6688 

8.0073% 

B  MFT-DHODA-  8 

0.2934 

0.0649 

0.3554 

0.1281% 

BMFT-DHODA- 16 

0.1677 

0.0412 

0.5608 

0.0813% 

HTDA 

0.5125 

0.0942 

0.3715 

0.0683% 

BMFT-HTDA-8 

0.2930 

0.0646 

0.1228 

0.0796% 

MSE(Oi,Di) 

MSE(D2,D2) 

RMSE 

DHODA 

0.6948 

12.9626 

1.83% 

B  MFT-DHODA-  8 

0.0427 

0.1310 

0.25% 

BMFT-DHODA- 16 

0.0419 

0.0701 

0.18% 

HTDA 

0.1276 

0.0430 

0.23% 

BMFT-HTDA-8 

0.0235 

0.0716 

0.20% 

Figure  1 1 :  Wideband  sinusoidal  FM  image  with  IF  components  outside  the  range  of  ESA  con¬ 
straint 


image  for  this  example  and  its  corresponding  IF  components  are  given  by: 

,2ti  2n 
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Note  that  the  carrier  frequency  along  each  direction  is  Hence  the  IF  components  of  the  im¬ 
age  are  out  of  the  range  [0,  |],  The  demodulation  results  for  the  DHODA  and  the  BMFT-DHODA 
are  illustrated  in  Fig.  12.  Specifically,  Fig.  12  contains  perspective  plots  of  the  IF  estimation  for  the 
wideband  sinusoidal  example  where  IF  components  are  out  of  the  range  [0,  |]  via  the  DHODA  and 
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y  x  y 


(C)  (d) 

Figure  12:  Perspective  plot  of  the  IF  estimation  for  the  wideband  sinusoidal  example  where  IF 
components  are  out  of  the  ESA  estimation  range  via  the  DHODA  and  via  the  BMFT-DHODA 


via  the  BMFT-DHODA  with  conversion  factors  R  =  diag[8,  8].  Fig.  12(a)  is  an  estimation  of  the 
IF  component  along  the  horizontal  direction  via  the  DHODA.  Fig.  12(b)  is  an  estimation  of  the  IF 
component  along  the  vertical  direction  via  the  DHODA.  Fig.  12(c)  estimation  of  the  IF  component 
along  the  horizontal  direction  via  the  BMFT-DHODA.  Fig.  12(d)  is  an  estimation  of  the  IF  com¬ 
ponent  along  the  vertical  direction  via  the  BMFT-DHODA.  As  we  can  observe,  the  demodulated 
IF  components  via  the  DHODA  are  seriously  distorted  in  both  amplitudes  and  phases,  whereas 
the  demodulated  IF  components  via  the  BMFT-DHODA  are  smooth  and  valid.  The  RMSE  for  the 
DHODA  is  as  large  as  56.74%  while  the  RMSE  for  the  BMFT-DHODA  is  merely  0.9113%.  Hence 
our  claim  that  the  BMFT  helps  overcome  the  range  constraint  of  demodulation  algorithms  based 
on  ESA  is  justified. 

In  the  third  experiment,  we  justify  the  efficacy  of  the  proposed  BMFT  approach  on  real  images. 
The  real  image  of  an  oakring  (photo  by  H.D.  Grissino-Mayer)  is  shown  by  Fig.  13(a)  where  (a)  is 
Real  image  of  the  oakring,  (b)  is  the  estimation  of  the  IF  component  along  the  horizontal  direction 
via  the  DHODA,  (c)  is  the  estimation  of  the  IF  component  along  the  vertical  direction  via  the 
DHODA,  (d)  is  the  estimation  of  the  IF  component  along  the  horizontal  direction  via  the  BMFT- 
DHODA,  (e)  is  the  estimation  of  the  IF  component  along  the  vertical  direction  via  the  BMFT- 
DHODA,  (f)  is  the  estimation  of  the  IF  needle  plot  via  the  DHODA,  (g)  is  the  estimation  of  the 
IF  needle  plot  via  the  BMFT-DHODA,  (h)  is  the  estimation  of  the  IA  via  the  DHODA,  and  (i)  is 
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(g)  (h)  (i) 

Figure  13:  Demodulation  of  the  oakring  image  (photo  by  H.D.  Grissino-Mayer) 


the  estimation  of  the  IA  via  the  BMFT-DHODA.  Note  that  the  conversion  factors  of  the  BMFT- 
DHODA  here  are  R  =  diag[16, 16].  The  demodulation  results  associated  with  this  real  oakring 
image  are  compared  in  Fig.  13.  The  estimated  IF  components  via  the  DHODA,  as  shown  in 
Fig.  13(b)  and  Fig.  13(c),  have  singular  points  with  significantly  large  values.  The  estimated  IF 
components  via  the  BMFT-DHODA,  as  shown  in  Fig.  13(d)  and  Fig.  13(e),  are  much  smoother 
and  do  not  have  any  singular  point  with  significantly  large  value.  As  a  result,  the  estimated  IF 
needle  plot  of  the  DHODA  (Fig.  13(f))  is  seriously  distorted  by  those  singular  points,  whereas  the 
estimated  IF  needle  plot  of  the  BMFT-DHODA  (Fig.  13(g))  reveals  the  ring  pattern  corresponding 
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to  the  real  oakring  image.  Note  that  the  IA  estimated  via  the  DHODA  (Fig.  13(h))  suffers  from 
the  same  issue  whereas  the  IA  estimated  via  the  BMFT-DHODA  (Fig.  13(i))  does  not.  We  see 
that  ESA  demodulation  employed  by  the  DHODA  prohibits  its  direct  application  to  wideband 
real  images  due  to  the  narrowband  and  estimation  range  constraints,  while  the  proposed  BMFT 
framework  allows  for  direct  application  of  ESA  demodulation  to  wideband  real  images. 
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4.2  Wideband  Hypernasal  Speech  Demodulation 

Formant  analysis  has  found  various  applications  in  fields  such  as  speech  synthesis,  bandwidth 
compression  and  speech  transmission.  Formants  are  natural  resonances  of  the  vocal  tract  and  each 
formant  is  generally  modeled  as  an  AM-FM  signal  [2].  The  goal  of  formant  analysis  here  is  to 
demodulate  the  IF  and  the  IA  of  the  corresponding  formant.  But  first  of  all,  we  need  to  isolate 
each  formant  from  the  speech  signal,  which  is  conventionally  achieved  via  bandpass  filtering. 
However,  by  using  bandpass  filtering,  the  center  frequency  and  bandwidth  are  manually  assigned 
to  each  formant,  which  may  incur  serious  distortion  to  the  isolated  resonance  signal  and  result  in 
inaccurate  tracking  of  the  IF  and  IA  of  the  corresponding  formant.  In  addition,  wideband  formants 
may  overlap  with  each  other  in  the  spectrum,  therefore  they  cannot  be  simply  isolated  via  bandpass 
filtering. 

In  this  section,  we  primarily  focus  on  the  first  formants  of  the  vowels,  which  are  generally 
wideband  due  to  the  small  center  (or  carrier)  frequency,  especially  the  nasal  vowels.  We  propose 
the  Empirical  Mode  Decomposition  (EMD)  [17]  instead  of  bandpass  filtering  for  the  first  formant 
extraction  and  apply  the  MFT  approach  to  demodulate  the  IA  and  IF  of  the  wideband  first  formant 
for  further  analysis. 

4.2.1  Empirical  Mode  Decomposition 

The  EMD  algorithm  is  an  intuitive  method  that  perform  decomposition  process  adaptively  with  an 
posteriori  defined  basis  derived  from  the  data  itself.  It  generally  involve  two  major  procedures, 
namely  the  sifting  process  and  decomposition. 

4.2.2  Sifting  Process 

A  function  is  called  an  intrinsic  mode  function  (IMF)  if  the  following  conditions  are  satisfied:  1) 
The  number  of  extremas  and  the  number  of  zero  crossings  equals  or  differs  at  most  by  one;  2)  the 
average  of  the  upper  envelope  defined  by  local  maximas  and  the  lower  envelope  defined  by  local 
minimas  at  any  point  is  zero.  The  IMF  reflects  the  oscillation  mode  inherent  in  the  signal.  It  can 
be  modeled  as  AM-FM  and  is  not  restricted  to  narrowband  signal. 

The  sifting  process  is  a  systematic  way  to  extract  the  IMF  from  the  input  data  x(t)  and  can  be 
summarized  via 

•  Initialize  d0(t)  =  x(t) 

•  Identify  the  local  extremas  of  dn(t). 

•  Interpolate  the  local  maximas  and  local  minimas  to  form  the  the  upper  envelope  un(t)  and 
lower  envelope  vn(t)  respectively. 

•  Determine  the  local  mean  of  the  upper  and  lower  envelopes  mn(t )  =  Un(t>+vn(t)  ^ 

•  Extract  the  detail:  dn+i(t )  =  dn(t )  —  mn(t). 

•  Repeat  from  step  2  to  step  5  until  dn+i(t)  is  an  IMF  (zero  mean  or  stopping  criterion  met). 
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4.2.3  Decomposition 


Assume  that  the  signal  X(t)  is  composed  of  a  few  oscillatory  modes  that  can  be  modeled  by  the 
IMF.  Decomposition  is  a  procedure  that  keep  repeating  the  sifting  process  to  decompose  the  signal 
as  the  sum  of  a  few  IMFs  plus  a  residue  described  by 


n 

X(t )  =  y>fc  +  rn,  (76) 

k=  1 

where  Ck(t)  denotes  the  corresponding  IMF  and  rn(t)  denote  the  final  residue. 

The  decomposition  procedure  is  summarized  via 

•  Initialize  ru(t)  =  X(t). 

•  Apply  the  sifting  process  on  rn(t)  to  obtain  the  corresponding  IMF  cn+1(f)  and  the  residue 

rn+i(t)  =  rn(t )  -  cn_|_i (t) . 

•  Repeat  the  previous  step  until  the  residue  rn+1  (t)  has  no  more  extremas  or  meets  the  stopping 
criterion. 

4.2.4  EMD  Analysis 

According  to  the  structure  of  the  EMD,  the  number  of  extrema  associated  with  the  extracted  IMF  is 
gradually  decreased  iterating  from  one  residue  to  the  next  in  the  decomposition  procedure.  It  corre¬ 
sponds  to  a  interbank  structure  with  subband  changing  from  high  frequency  range  to  low  frequency 
range.  However,  it  is  different  from  any  predetermined  subband  filtering,  since  the  frequency  range 
and  resolution  associated  with  each  subband  are  adaptively  time- varying.  It  offers  more  flexibility 
than  the  traditional  multiband  filtering  approach  in  capturing  features  that  are  nonstationary. 

Except  the  basic  EMD  algorithm  reviewed  in  previous  sections,  there  are  a  few  variants  such  as 
the  ensemble  empirical  mode  decomposition  (EEMD)  [19].  The  key  idea  of  the  EEMD  relies  on 
averaging  the  modes  obtained  by  applying  EMD  to  different  realizations  of  Gaussian  white  noise 
added  to  the  original  signal.  By  introducing  white  Gaussian  noise,  the  mode  mixing  problem  can 
be  solved  via  populating  the  whole  time-frequency  space  to  take  advantage  of  the  dyadic  filter  bank 
behavior  of  the  EMD. 

Currently,  we  only  employ  the  basic  EMD  to  obtain  the  first  formant  of  the  vowel  for  the 
purpose  of  illustration.  In  the  future,  we  intend  to  apply  other  variants  of  the  EMD  algorithm  such 
as  the  EEMD  with  adaptive  noise  to  see  if  better  results  can  be  achieved. 

4.2.5  First  Formant  Extraction  via  EMD 

Formants  are  natural  resonances  of  the  vocal  tract  and  have  considerable  amounts  of  amplitude  and 
frequency  modulation.  Speech  quality  may  deteriorate  significantly  when  the  modulations  are  re¬ 
moved.  The  center  frequency  and  bandwidth  of  the  formant  associated  with  any  given  vowel  varies 
from  person  to  person.  In  addition,  the  first  formant  usually  lies  in  the  low-frequency  range  and 
has  a  relatively  small  center  frequency  to  bandwidth  ratio.  Therefore  it  is  generally  considered  as  a 
wideband  formant.  In  this  effort,  the  EMD  is  preferred  over  bandpass  filtering  as  the  technique  to 
extract  the  first  formant  due  to  the  following  reasons:  1)  EMD  does  not  require  accurate  estimation 
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for  the  center  frequency  and  bandwidth  of  the  first  formant,  which  is  hardly  accessible  via  Linear 
Predictive  Coding  (LPC)  or  other  methods.  2)  EMD  allows  for  some  amount  of  spectrum  overlap 
between  the  first  formant  and  the  second  formant. 

To  extract  the  first  formant,  we  perform  the  EMD  on  the  speech  signal  S(t)  to  obtain  the  IMFs 
and  the  residue  through 

n 

S(t )  =  y^Cfc  (t )  +  r  n  ( t ) .  (77) 

k= 1 

By  choosing  an  appropriate  number  m,  we  can  ignore  the  first  m  IMFs  that  locate  in  the  high- 
frequency  range  of  the  spectrum  and  sum  up  the  rest  of  the  IMFs  plus  the  residue  to  obtain  the  first 
formant  via 

n 

F\{t)=  ^2  Ck(t)  +  rn(t).  (78) 

k=m-\- 1 

We  can  see  that  the  first  formant  extraction  via  the  EMD  does  not  require  any  estimation  for  the 
center  frequency  and  the  bandwidth.  It  is  simple  to  choose  an  appropriate  number  m  by  observing 
the  oscillation  mode  of  each  IMF.  In  addition,  since  each  IMF  is  modeled  as  AM-FM  signal,  some 
amount  of  spectrum  overlapping  between  the  first  formant  and  the  second  formant  is  allowed. 
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4.2.6  Example  of  The  First  Formant  Demodulation  via  the  MFT  approach 


Figure  14:  Example  of  the  wideband  first  formant  extraction  via  EMD 


In  this  section,  we  present  an  example  of  the  demodulation  of  the  first  wideband  formant.  The 
speech  waveform  and  the  corresponding  short-time  spectrum  of  an  women’s  nasal  vowel  /iy/  (sam¬ 
ple  collected  by  Christopher)  is  shown  by  Fig.  14(a)  and  Fig.  14(b)  respectively.  The  first  formant 
of  this  vowel  is  obviously  wideband,  with  a  relatively  small  center  frequency  but  a  large  bandwidth. 
The  short-time  spectrum  of  the  first  formant  extracted  via  the  EMD  as  we  previously  propose  is 
shown  by  Fig.  14(c).  For  comparison,  the  short-time  spectrum  of  the  first  formant  extracted  via 
bandpass  filtering  using  FPC  estimation  of  center  frequency  and  bandwidth  is  shown  by  Fig.  14(d). 
Since  the  FPC  estimation  approach  is  only  valid  for  narrowband  signals,  the  estimated  bandwidth 
of  the  first  formant  is  very  small,  thus  leading  to  inaccurate  extraction  of  the  first  formant.  On  the 
other  hand,  the  distortion  of  the  first  formant  extracted  via  the  EMD  is  very  small  compared  to  the 
original  spectrum. 
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Instantaneous  Frequency 


Figure  15:  IF  demodulation  of  the  first  formant 
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Since  the  first  formant  in  this  example  is  wideband,  we  apply  the  MFT  and  Hilbert  transform 
combination  as  the  demodulation  scheme  to  obtain  the  IF.  The  IF  demodulation  of  the  first  formant 
is  illustrated  in  Fig.  15.  The  blue  line  corresponds  to  the  IF  demodulated  via  MFT  and  Hilbert 
transform  using  the  first  formant  extracted  via  the  EMD.  For  comparison,  the  red  line  corresponds 
to  the  IF  demodulated  via  the  Hilbert  transform  alone  using  the  first  formant  extracted  via  band¬ 
pass  filtering  in  conjunction  with  LPC  estimation  of  the  center  frequency  and  bandwidth.  As  we 
can  observe,  the  IF  demodulated  from  the  first  formant  extracted  via  bandpass  filtering  is  almost 
a  straight  line,  reflecting  the  fact  that  the  bandpass  filtering  with  LPC  estimation  only  works  un¬ 
der  the  narrowband  scenario,  whereas  the  IF  demodulated  using  the  MFT  and  Hilbert  transform 
combination  keeps  allows  for  tracking  the  IF  variation  of  the  first  formant  in  a  much  wider  range. 

4.3  Wideband  CPM  Demodulation  for  Satellite  Communications 

Wideband  FM  demodulation  is  a  popular  technique  for  satellite  communications  due  to  its  ability 
to  deal  with  trans-ionospheric  distortion  [3].  Contemporary  methods  for  the  demodulation  of  wide¬ 
band  FM  signals  either  use  negative  feedback  [2]  or  frequency  estimation,  or  multirate  frequency 
transformations  and  adaptive  frequency  tracking  [4] .  A  particular  form  of  digital  FM,  multilevel 
Gaussian  Frequency  Shift  Keying  (FSK)  has  been  proposed  as  a  solution  for  high  bandwidth  satel¬ 
lite  communications. 

On  the  wireless  communication  front,  continuous  phase  modulation  (CPM)  belongs  to  a  class 
of  non  linearly  modulated  signals  with  constant  envelope,  where  the  information  is  carried  in 
the  phase  of  the  transmitted  signal.  High  spectral  efficiency  and  suitability  to  non  linear  class  C 
amplifiers  used  in  mobile  radio  applications  make  CPM  a  popular  modulation  choice.  A  specific 
form  of  CPM  namely  Gaussian  minimum  shift  keying  (GMSK)  has  been  adopted  in  the  Global 
System  for  Mobile  communications  (GSM)  [21,  22].  The  optimum  receiver  structure  for  CPM 
demodulation  employs  the  maximum  likelihood  (ML)  detector  based  on  the  Viterbi  algorithm 
[18,  21,  22].  This  receiver  structure,  however  has  significant  computational  complexity  which 
grows  exponentially  with  increase  in  the  number  of  phase  states. 

In  prior  work  [5,  8],  it  was  shown  that  frequency  discrimination  for  full  response  CPM  demod¬ 
ulation  has  the  same  performance  as  that  of  binary  phase  shift  keying  (BPSK)  detection  in  AWGN. 
In  recent  work  [4]  frequency  tracking  based  wideband  FM  demodulation  was  extended  to  large 
wideband  to  narrowband  conversion  factors  using  multirate  frequency  transformations.  Frequency 
estimation  based  approaches  have  the  added  advantage  that  they  are  immune  to  phase  distortions 
introduced  by  the  channel  which  would  affect  the  Vitterbi  approach.  In  this  section,  we  apply  the 
MFT  enabled  wideband  FM  demodulation  approach  to  CPM  demodulation.  The  first  step  is  the 
adaptive  linear  predictive  tracking  of  the  instantaneous  frequency  of  the  CPM  waveform  using  a 
gradient  based  adaptive  step-size  normalized  least  mean  squares  (LMS)  algorithm  [7]  and  multirate 
frequency  transformations. 

4.3.1  CPM  Signal  Model 

For  discussion  purposes,  we  will  adopt  a  rectangular  pulse-shaping  function  p(t)  with  a  duration 
of  L  symbol  periods  and  binary  Pulse  amplitude  Modulation  (PAM)  symbols  a[k]  G  {  —  1,1}.  The 
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IF  demodulation  of  the  CPFSK  waveform 


LPC  Error 


(c) 


Figure  16:  Two-stage  linear  prediction  approach  for  CPFSK  demodulation 


IF  signal  in  either  case  takes  the  form: 

OO 

0Ji(t)  =  loc  +  2nh  a[k]p(t  —  kTb) ,  (79) 

k=—oo 

where  uic  is  the  carrier  frequency  and  h  is  the  modulation  index  of  CPM.  The  phase  deviation  from 
the  carrier  phase  is  given  by: 


OO 

4>dev{t\  a)  =  2,7th  ^  a[k\q(t  -  kTb ), 

k=— oo 


(80) 


where  q(t)  =  f*p(r)dT  corresponds  to  the  phase  pulse  shaping  function.  The  CPM  signal  is  then 
obtained  via  frequency  modulation: 


r(t)  =  A  cos 


u}i(r)dT  +  60 


(81) 
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Using  a  pulse  shaping  function  of  duration  larger  than  a  symbol  period  introduces  memory  into  the 
modulation  scheme  (LREC-CPM).  In  this  report,  we  will  focus  our  attention  on  the  memory  less 
case  with  L  —  1,  i.e.,  (1REC-CPM).  Specifically  CPM  with  a  rectangular  pulse  of  one  symbol  du¬ 
ration  (1-REC-CPM)  is  equivalent  to  continuous  phase  frequency  shift  keying  (CPFSK).  Another 
form  of  CPM,  minimum  shift  keying  (MSK),  is  equivalent  to  1-REC-CPM  with  h  =  0.5,  while 
GMSK  can  be  put  into  the  CPM  framework  with  a  Gaussian  pulse  function.  The  results  of  which 
are  shown  in  Fig.  16  where  (a)  is  the  modulating  bits  and  corresponding  CP-FSK  waveform  for 
2-REC  CPM  at  a  SNR  of  13  dB,  (b)  is  the  IF  estimate  from  the  first  stage  of  adaptive  linear  predic¬ 
tion,  and  (c)  is  the  removal  of  memory  induced  by  partial  response  signaling  using  a  second-stage 
of  linear  prediction. 
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4.3.2  Wideband  CPM  Demodulation 


f  =  10  kHz,  T  =  0.5  ms,  n  =  500,  h  =  0.1,  REC-CPM 
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Figure  17:  Performance  of  the  dual  stage  linear  prediction  approach 


In  prior  work  [4],  we  have  applied  various  forms  of  adaptive  filtering  to  the  demodulation  of 
CPM  signals  without  memory.  It  was  shown  that  oversampling  results  in  an  improved  performance 
analogous  to  that  of  binary  detection  in  AWGN.  When  partial  response  signaling  is  employed,  this 
induces  ISI  in  the  estimated  IF  and  a  loss  of  performance.  Figure  17  depicts  the  application  of 
the  wideband  demodulation  approach  to  CPM  signals  with  memory,  i.e.,  L  =  2  CPM.  As  can  be 
observed,  partial  response  signaling  induces  ISI  and  a  deterioration  in  performance  after  averaging 
over  500  experiments.  Specifically  adjusting  the  predictor  order  allows  for  removal  of  more  of  the 
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memory  induced  ISI.  Within  Fig.  17,  (a)  is  for  L  =  2-REC  CPM  using  a  p  =  20-th  order  predictor 
Partial  response  signaling  induced  ISI  results  in  a  error  floor  that  can  be  mitigated  to  a  certain  extent 
by  the  second  stage  of  linear  prediction,  (b)  is  the  effect  of  predictor  order  in  the  second  stage  on 
the  performance  of  the  algorithm,  and  (c)  is  the  effect  of  pulse  shape  with  no  memory.  Specifically 
future  work  will  focus  on  the  following:  (a)  the  effect  of  channel  memory  on  the  performance  of 
the  proposed  approach,  (b)  effect  of  predictor  order  on  performance,  and  (c)  effect  of  pulse  shape 
on  the  performance  of  the  algorithm.  Equalization  efforts  such  as  decision  feedback  equalization 
[18]  applied  directly  to  the  estimated  IF  can  be  used  to  further  aid  in  the  reduction  of  the  ISI. 
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5  Conclusion 


In  this  report,  we  have  presented  an  extension  of  the  recently  proposed  [4]  MFT  based  wideband 
FM  demodulation  approach  to  large  wideband  to  narrowband  demodulation  factors,  thereby  en¬ 
abling  more  practical  implementations  of  the  filters  internal  to  the  heterodyne  and  interpolation 
modules.  It  was  shown  that  this  approach  significantly  reduces  the  demodulation  errors  in  compar¬ 
ison  to  direct  application  of  the  narrowband  approaches. 

We  further  extended  the  MFT  approach  to  two  dimensions  and  combined  it  with  demodulation 
approaches  based  on  higher-order  two-dimensional  Teager- Kaiser  energy  operators.  When  applied 
to  both  synthetic  and  natural  wideband  images  to  demonstrate  that  significant  reduction  of  demod¬ 
ulation  errors  over  just  narrowband  demodulation  algorithms.  When  applied  to  wideband  formants 
encountered  in  hypernasal  situations,  the  wideband  demodulation  approach  provides  more  accu¬ 
rate  estimates  of  the  instantaneous  frequency  of  the  first  formant  over  linear  prediction  approaches 
that  filter  out  a  lot  of  the  valuable  sidelobe  information. 

When  then  applied  the  proposed  approach  to  the  problem  of  wideband  CPM  demodulation, 
the  approach  with  oversampling  was  shown  to  provide  a  performance  similar  to  that  of  binary 
detection  in  AWGN  in  the  memoryless  situation.  For  symbols  with  memory,  linear  prediction 
of  the  IF  estimates  removes  some  of  the  induced  memory,  but  further  equalization  in  the  form  of 
minimum  mean  squared  error  or  decision  feedback  equalization  [1 8]  of  the  estimated  instantaneous 
frequency  is  needed  to  completely  remove  the  memory. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS  AND  ACRONYMS 


AF-RLS 

Root  Feast  Squares 

AM-FM 

Amplitude  modulation  frequency  modulation 

AS 

Analytic  signal 

AS-LMS 

Analytic  Signal  Feast  Mean  Square 

AWGN 

Additive  White  Gaussian  Noise 

BEMD 

Bi-dimensional  Empirical  Mode  Decomposition 

BMFT 

Bi-dimensional  multirate  frequency  transformation 

BPF 

Band  Pass  Filter 

BPSK 

Binary  Phase  Shift  Keying 

CPFSK 

Continuous-Phase  Frequency-Shift  Keying 

CPM 

Continuous  phase  modulation 

CR/FD 

Carrier  to  frequency  deviation  ratio 

CR/IB 

Carrier  to  information  bandwidth  ratio 

CRLB 

Cramer-Rao  lower  bound 

dB 

Decibel 

DCA 

Dominant  Component  Analysis 

DHODA 

Discrete  higher  order  demodulation  algorithm 

EEMD 

ensemble  empirical  mode  decomposition 

EMD 

Empirical  Mode  Decomposition 

ESA 

Energy  separation  algorithm 

FIR 

Finite  Impulse  Response 

FSK 

Frequency  Shift  Keying 

GMSK 

Gaussian  minimum  shift  keying 

GNGD 

Generalized  normalized  gradient  descent 

HT 

Hilbert  transform 

HTDA 

partial  Hilbert  transform  alone 

IA 

Instantaneous  amplitude 

IF 

Instantaneous  frequency 

IFGT 

instantaneous  frequency  gradient  tensor 

IMF 

intrinsic  mode  function 

FMS 

Feast  Mean  Squares 

FPC 

Finear  predictive  coding 

MFT 

Multirate  frequency  transformation 

MSK 

Minimum  shift  keying 

NFMS 

Normalized  Feast  Mean  Squares 

NRMSE 

Normalized  Root  Mean  square  error 

PAM 

Pulse  Amplitude  Modulation 

RMSE 

Root  Mean  Square  Error 

SNR 

Signal  to  Noise  Ratio 
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